/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://github.com/shenyczz/CSharpKit
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/
using System;

namespace CSharpKit.SVM
{
    /// <summary>
    /// 核
    /// </summary>
    internal abstract class Kernel : IQMatrix
    {
        /// <summary>
        /// 构造函数
        /// </summary>protected
        /// <param name="l">样本集长度</param>
        /// <param name="x">样本节点</param>
        /// <param name="param">核函数参数</param>
        protected Kernel(int l, SvmNode[][] x, SvmParameter param)
        {
            _kernelType = param.KernelType;
            _degree = param.Degree;
            _gamma = param.Gamma;
            _coef0 = param.Coef0;

            _x = (SvmNode[][])x.Clone();

            if (_kernelType == KernelType.RBF)
            {
                _xSquare = new double[l];
                for (int i = 0; i < l; i++)
                    _xSquare[i] = dot(_x[i], _x[i]);
            }
            else
            {
                _xSquare = null;
            }
        }


        #region Fields

        private KernelType _kernelType;
        private int _degree;
        private double _gamma;
        private double _coef0;

        private SvmNode[][] _x;
        private double[] _xSquare;  // ?

        #endregion


        #region IQMatrix Members

        public abstract float[] GetQ(int i, int len);
        public abstract double[] GetQD();
        public virtual void SwapIndex(int i, int j)
        {
            _x.SwapIndex(i, j);

            if (_xSquare != null)
            {
                _xSquare.SwapIndex(i, j);
            }
        }

        #endregion


        /// <summary>
        /// 核函数
        /// </summary>
        /// <param name="i"></param>
        /// <param name="j"></param>
        /// <returns></returns>
        public double KernelFunction(int i, int j)
        {
            double result = 0;

            switch (_kernelType)
            {
                case KernelType.LINEAR:
                    result = dot(_x[i], _x[j]);
                    break;
                case KernelType.POLY:
                    result = powi(_gamma * dot(_x[i], _x[j]) + _coef0, _degree);
                    break;
                case KernelType.RBF:
                    result = System.Math.Exp(-_gamma * (_xSquare[i] + _xSquare[j] - 2 * dot(_x[i], _x[j])));
                    break;
                case KernelType.SIGMOID:
                    result = Math.Tanh(_gamma * dot(_x[i], _x[j]) + _coef0);
                    break;
                case KernelType.PRECOMPUTED:
                    result = _x[i][(int)(_x[j][0].Value)].Value;
                    break;
                default:
                    result = 0;
                    break;
            }

            return result;
        }

        public static double KernelFunction(SvmNode[] x, SvmNode[] y, SvmParameter param)
        {
            double result = 0;

            switch (param.KernelType)
            {
                case KernelType.LINEAR:
                    result = dot(x, y);
                    break;
                case KernelType.POLY:
                    result = powi(param.Degree * dot(x, y) + param.Coef0, param.Degree);
                    break;
                case KernelType.RBF:
                    double sum = computeSquaredDistance(x, y);
                    result = System.Math.Exp(-param.Gamma * sum);
                    break;

                case KernelType.SIGMOID:
                    result = Math.Tanh(param.Gamma * dot(x, y) + param.Coef0);
                    break;
                case KernelType.PRECOMPUTED:
                    result = x[(int)(y[0].Value)].Value;
                    break;
                default:
                    result = 0;
                    break;
            }

            return result;
        }


        /// <summary>
        /// 点积
        /// </summary>
        /// <param name="xNodes"></param>
        /// <param name="yNodes"></param>
        /// <returns></returns>
        private static double dot(SvmNode[] xNodes, SvmNode[] yNodes)
        {
            double sum = 0;
            int xlen = xNodes.Length;
            int ylen = yNodes.Length;
            int i = 0;
            int j = 0;
            SvmNode x = xNodes[0];
            SvmNode y = yNodes[0];
            while (true)
            {
                if (x.Index == y.Index)
                {
                    sum += x.Value * y.Value;
                    i++;
                    j++;
                    if (i < xlen && j < ylen)
                    {
                        x = xNodes[i];
                        y = yNodes[j];
                    }
                    else if (i < xlen)
                    {
                        x = xNodes[i];
                        break;
                    }
                    else if (j < ylen)
                    {
                        y = yNodes[j];
                        break;
                    }
                    else break;
                }
                else
                {
                    if (x.Index > y.Index)
                    {
                        ++j;
                        if (j < ylen)
                            y = yNodes[j];
                        else break;
                    }
                    else
                    {
                        ++i;
                        if (i < xlen)
                            x = xNodes[i];
                        else break;
                    }
                }
            }
            return sum;
        }

        /// <summary>
        /// 浮点数的整数次幂
        /// </summary>
        /// <param name="value">底数</param>
        /// <param name="times">幂指数</param>
        /// <returns></returns>
        private static double powi(double value, int times)
        {
            double tmp = value, ret = 1.0;

            for (int t = times; t > 0; t /= 2)
            {
                if (t % 2 == 1) ret *= tmp;
                tmp = tmp * tmp;
            }

            return ret;
        }

        private static double computeSquaredDistance(SvmNode[] xNodes, SvmNode[] yNodes)
        {
            SvmNode x = xNodes[0];
            SvmNode y = yNodes[0];

            int xLength = xNodes.Length;
            int yLength = yNodes.Length;
            int xIndex = 0;
            int yIndex = 0;
            double sum = 0;

            while (true)
            {
                if (x.Index == y.Index)
                {
                    double d = x.Value - y.Value;
                    sum += d * d;
                    xIndex++;
                    yIndex++;
                    if (xIndex < xLength && yIndex < yLength)
                    {
                        x = xNodes[xIndex];
                        y = yNodes[yIndex];
                    }
                    else if (xIndex < xLength)
                    {
                        x = xNodes[xIndex];
                        break;
                    }
                    else if (yIndex < yLength)
                    {
                        y = yNodes[yIndex];
                        break;
                    }
                    else break;
                }
                else if (x.Index > y.Index)
                {
                    sum += y.Value * y.Value;
                    if (++yIndex < yLength)
                        y = yNodes[yIndex];
                    else break;
                }
                else
                {
                    sum += x.Value * x.Value;
                    if (++xIndex < xLength)
                        x = xNodes[xIndex];
                    else break;
                }
            }

            for (; xIndex < xLength; xIndex++)
            {
                double d = xNodes[xIndex].Value;
                sum += d * d;
            }

            for (; yIndex < yLength; yIndex++)
            {
                double d = yNodes[yIndex].Value;
                sum += d * d;
            }

            return sum;
        }




        //}}@@@
    }



}
